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Field Of The Invention 

The present invention relates generally to imaging, and more particularly to 
a physically-based, probabilistic model for ultrasonic images incorporating shape, 
microstructure and system characteristics. 

Background Of The Invention 

Necessary in the art of image modeling is an accurate, physically-based, 
computationally feasible probabilistic model for ultrasonic images. Model-based image 
analysis has become increasingly popular and successful to the analysis of object shape in 
images. Traditional approaches, usually termed pattern or object recognition, typically 
involve building a pattern or object from features extracted by "processing" the image. In 
the model-based approach, shape is modeled rigorously then inferred using an image 
model. Shape can be inferred in a Bayesian setting by combining a probabilistic prior 
model describing variation of the shape with a data likelihood (i.e., a probabilistic model 
describing observations of the shape). The physics of the imaging system are incorporated 



using the data likelihood. Model-based approaches can provide insight and understanding 
regarding the process of image formation and its connection to underlying structure, 
thereby forming a solid foundation for inference of structural patterns in image data. 

Shape is represented in ultrasound as the result of complex interactions at the 
5 microstructural level, producing a mix of speckle texture and coherent echoes. Application 
of model-based techniques, therefore, require a probabilistic model unique to ultrasound 
that accurately represents shape in terms of these interactions. Previous ultrasound 
applications of model-based approaches to shape analysis have been rmnimal and have used 
: simplistic data models (e.g., Rayleigh with a constant parameter over a given contour). 
10 t- Previous applications have shown some success with simple shape models, although only in 
* high contrast images. Because previous applications are not physically-based, the results 
are not extendable to other image analysis problems and little intuition or insight is 
developed for representing and understanding the relation between the data and the 
l- underlying shape. 

15 For the foregoing reasons, there is a need for an accurate, physically-based, 

computationally-feasible probabilistic model for ultrasonic images incorporating shape, 
microstructure and system characteristics. 

Summary of the Invention 
The present invention is directed to a method and apparatus for forming an 
20 image model and for developing a physically-based, probabilistic model for ultrasonic 
images incorporating shape, microstructure and system characteristics. The present 
invention accurately represents complex shape, allowing inference of underlying structural 
patterns in the image data. 
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An aspect of the present invention provides a method and an apparatus to 
perform the method for developing a physical model of image formation and performing a 
random phasor sum representation of the physical model to form a probabilistic model. 

In one aspect of the present invention, forming the physical model includes 
5 formulating imaging system characteristics, formulating gross shape, formulating 
microstructure and incorporating the imaging system characteristics, the gross shape and 
the microstructure to form the image model. 

In another aspect of the present invention, the imaging system characteristics 
can be formulated using a three-dimensional point spread function. The imaging system 
10 characteristics could also include a radio frequency image described by a linear model and 
characterized by the point spread function. 

In yet another aspect of the present invention, the image model includes a 
data likelihood enabling a statistical inference to formulate underlying characteristics. The 
data likelihood can be formulated using image pixel based statistics. 
15 In still another aspect of the present invention, using the image pixel based 

statistics can include computing an amplitude mean value, computing an amplitude variance 
value and computing a ratio of the amplitude mean value to a standard deviation value at 
each image pixel to develop a statistical image characterizing tissue. Then, classifying each 
pixel as Rayleigh or non-Rayleigh determined by the ratio of the amplitude mean value to 
20 the standard deviation value. Assigning a density function to each image pixel based upon 
the classification of each image pixel and then constructing the data likelihood as a product 
of the density functions. 

In a further aspect of the present invention, the gross shape is described by a 
triangulated surface which can include a set of triangular elements defined by respective 
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vertices and edges. Acoustic properties of the triangulated surface could be represented by 
multiple discrete scatterers distributed across the triangulated surface in a random model. 
Spatial locations of the scatterers across the triangulated surface could be parametrized by a 
scatterer concentration and a surface roughness. 
5 In a still further aspect of the present invention, each discrete scatterer could 

contribute a phasor to the random phasor sum representation of the physical model. 

It is therefore an object of the present invention to develop a probabilistic 
image model for ultrasound applicable to the Bayesian, model-based approach to image 
•'■ analysis. 

10 r: Another object of the present invention is to include the gross shape, the 

surface microstructure and the imaging system point spread function into the probabilistic 
model for ultrasound. 

.:■ Still another object of the present invention is to develop a probabilistic 

model for ultrasound requiring minimal computational demands while providing accurate 
15 characterization of the pixel intensity variation across the image. 

Yet another object of the present invention is to develop a probabilistic 
model for ultrasound providing a representation suitable for construction of a data 
likelihood for the image offering smooth variation for changes in the object shape making 
the probabilistic model suitable to derivative based optimization techniques for inference. 
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Brief Description Of The Drawings 

For the purpose of illustrating the invention, there is shown in the drawings a 
form which is presently preferred; it being understood, however, that this invention is not 
limited to the precise arrangements and instrumentalities shown. 
5 Figure 1 illustrates a phantom containing a lumbar vertebra constructed to 

allow registration of CT images, a vertebra triangulated surface representation and 
ultrasonic images in accordance with the present invention; 

Figure 2 illustrates approximate planes for the development of images in 
accordance with the present invention, the planes overlayed on a rendering of a triangulated 
ion model for the vertebral surface with transverse process shown on the left and lamina shown 
on the right; 

JL, Figure 3 illustrates an actual image at top left and three simulated images at 

y bottom left of the transverse process with corresponding statistical images at right, the 
CI images developed in accordance with the present invention; 
15 Figure 4 illustrates an actual image at the top and three simulated images 

below of the lamina and articular processes, the images developed in accordance with the 

present invention; 

Figure 5 illustrates statistical images of the lamina and articular processes 
with regions of non-Rayleigh scattering existing along the lamina and peak of the articular 
20 process, the images developed in accordance with the present invention; and 

Figure 6 illustrates samples of a Rayleigh/Guassian image model of the 
lamina with pixel intensities generated using computations in accordance with the present 
invention. 
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Detailed Description Of The Invention 

A probabilistic image model of the present invention is derived from a phys- 
ical model for image formation as described in section I. The physical model describes the 
image formation process given a deterministic description of the gross shape, a random 
5 description of the associated microstructure, and a deterministic description of the system 
characteristics. For any realization of the random microstructure, a simulated image can be 
generated describing a "typical" system output. The probabilistic model of the present 
invention is based on a random phasor sum representation of the physical model and char- 
=' acterizes the ensemble of images that can be generated from the physical model. Statistics 
10 l" (mean and variance) are computed directly for each pixel and then used to assign a density 
- function to each pixel for construction of a data likelihood (i.e., a probabilistic model 
describing observations of shape). Computation of the mean and variance is nontrivial, 
f requiring several surface integrals for each pixel. General solutions for the mean and 
= variance are developed in section II. Solutions specific to surfaces, including numerical 
15 approximations, are shown in section III. 

The description of the present invention focuses on surfaces because, as a 
descriptor of shape, surfaces represent an important aspect of tissue with a range of 
ultrasonic characteristics presenting a significant challenge in image modeling. Statistics 
computed in the model are shown for a cadaveric vertebra. The vertebral surface provides 
20 a good medium for evaluation because of its intricate curvature and sub-wavelength 
roughness. Imaging of the vertebrae is also of interest in the area of treatment guidance 
(delivery of radiosurgery and guidance of traditional surgery) based on CT images of the 
spine. Ultrasonic images can provide an estimate of position and orientation of a vertebra 
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relative to pre-operative CT images. This approach can lead to less invasive surgery and 
radiosurgery. 

I. A PHYSICAL MODEL OF IMAGE FORMATION 
The physical model for image formation provides a basis for the probabilistic 
5 image model of the present invention. A radio frequency (RF) portion of the imaging 
system is described by a linear model with an imaging system characterized by its point 
spread function (PSF). Gross surface shape is described in detail with a triangulated 
surface, as shown in Figure 2. Acoustic properties of the surface are represented entirely 
~ by a collection of discrete scatterers. A random model describes spatial locations of the 
10 I scatterers across a spatial extent of the surface (position is assumed to be uniformly distrib- 
- uted on the surface), and the surface is parametrized by a scatterer concentration 

(scatterers/mm 2 ) and a surface roughness. 
; A. Imaging System Model 
15 :: A 3-dimensional (3D) volume of image data, irf (x, y, z), is modeled as an 

output of a linear system with 3D impulse response, or point-spread function (PSF), h(x, y, z), 
to an input of the 3D reflectivity function, r(x, y, z). The reflectivity function is assumed to 
represent a local change in density and compressibility. Note that a density term can be 
included in the reflectivity function by assuming that only direct reflection (backscatter) is 
20 measured. 

irf (x, y,z)=h (x, y,z)*r (x, y, z) (1 ) 

An exact description of the PSF, which characterizes the system, requires 
numerical modeling based on geometry of the transducer, but simpler descriptions are used to 
characterize the system when the exact PSF is unknown or unnecessary. A typical 
25 simplifying choice for the PSF h(x, y, z) is a 3D Gaussian envelope modulating a sinusoidal 
waveform in the axial, z, direction, of wavenumber, k 0 , e.g., 
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h(x, y, z) = A(x, y, z) cos(2k 0 z) (2) 

A(x,y,z,) = e ff ' 2 e° y e a * (3) 
5 where the standard deviations, a x , o y ,a z , denote the beam widths in each of the directions. 

For comparison to actual images, the envelopes of the RF images were 
generated along the axial dimension. Envelope detection was accomplished by taking the 
magnitude of the complex signal, generated with the Hilbert transform. 
B. Surface Model 

10 . B The surface was modeled as a collection of discrete scatterers, a representation 

"j: having many advantages over a continuous representation. These advantages stem from a 
::: simplified view of the combined effects of system and tissue characteristics. In the discrete - 
\ scatter representation, each scatterer represents a major scattering element of size smaller 
than a wavelength. Roughness of surfaces is typically characterized in such a way that small 
15 (sub-wavelength) hills and valleys cover the surface. A hill or valley pointing in the direction 
~ of the transducer could be considered a major scattering element in the model of the present 
invention. A collection of small scattering elements comprises the entire acoustic 
representation of the surface. 

Discrete representation offers the following advantages: 
20 • Computation is simplified because the small scattering elements enable a unification of 
structure and imaging system in a linear systems approach. For a continuous representation, 
the surface integral equation for the scattered wavefield requires numerical computation for 
each frequency of the incident wavefield and summation to calculate the image. 
• A similar intuitive simplicity results from considering the surface as a collection of 
25 distinct elements, with their contributions adding constructively or destructively depending 
on phase separation, rather than considering the solution to a surface integral equation. 



-8- 



• The discrete-scatterer model allows use of the convenient random-walk analysis of 
scattering. Most probabilistic models are based on the random-walk analysis model, 
providing an extensive background of previous work for the development of a new 
probabilistic image model. 

5 The gross surface is represented as a triangulated surface (i.e., a set of 

triangular elements defined by their vertices and edges). The triangulated surface 
representation can be easily produced from a volume segmentation using the Marching Cubes 
algorithm. While only the triangles are used, the triangulated surface representation 
inherently defines second-order (curvature) information, which can be used to accurately 

10 describe the surface. Several computational algorithms exist in the computer graphics arena 
for the manipulation of these surfaces. 

The discrete version of the triangulated surface permits a wide variety of 
- medium parametrization. Scatterer spacing on the surface can be characterized in many 
ways, from completely regular spacing to random spacing. In addition to spacing, the surface 

15 _ can be characterized by scatterer concentration. Surface roughness can also be incorporated 
independently as a perturbation of each scatterer in a direction normal to the surface. 

Any choice for distributing scatterers on the surface results in a discrete 
representation consisting of a collection X = {(x h y„ z$, i e [1,.., N]} of N scatterers, with 
scatterer i of amplitude A l at position (x„ y„ z f ). The reflectivity function consists of a sum of 

20 appropriately scaled 3D delta functions, 8(x h y i? zj, 

r{x, y. z) = ]T A5(j: - x, . y - y u z - z,) (4) 

The RF image is then a sum of a scaled and a delayed version of the PSF, 

i rf (x.y,z) = h(x,y,z)*r(x,y,z) (5) 
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In the frequency domain, the convolution is a product allowing the tissue response to be 
represented as a sum of scaled delay terms, 

Irf(u,v,w) = H(u.v,w)R(u,v,w) (7) 
= H(u,v,tv)jr A >e- y[ul * +vy ' +wz ' ] (8) 

5 

The system response is assumed bandlimited, thus the image response can be 
computed exactly for any number of scatterers, without the limitation on scatterer positions 
imposed by a uniform grid approach. Computation can be performed in the spatial domain as 

10 a sum of shifted versions of the PSF. Alternatively, a complete tissue response can be 
computed in the frequency domain with the image computed using an inverse transform of 
the product of tissue and system responses. 

Inherent in the model are a number of assumptions. In the linear model, the 
Born approximation is assumed. It is also assumed that a similar representation would be 

15 valid for other gross shapes (e.g., volumes, curves and points), which are necessary for 
describing a typical tissue medium. For the vertebra, the bony surface is assumed to com- 
pletely attenuate the beam, thus the surface is currently processed to remove occluded 
elements from view. Soft tissue surfaces are expected to be similar in basic appearance 
(i.e., presence of coherent echoes and speckle texture). Note that for a tissue comprised of 

20 a surface surrounding a volume of scattering structure, an additional model would be 
included to represent the interior volume. 

II. AMPLITUDE MEAN AND VARIANCE VIA THE 
RANDOM PHASOR SUM 

25 The image model requires a probability density function describing the echo 

amplitude at each image pixel. In the present invention, the amplitude is characterized by 

its mean and variance at every image pixel with comprehensive treatment of the shape, 
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microstructure and system characteristics. The derivation begins with a random phasor 
sum representation of the physical model of image formation. 

A. Image Formation as a Random Phasor Sum Representation of the Physical Model 

The phasor sum model of the present invention results from a linear model 
5 for a radio frequency (RF) image, where the imaging system is characterized by its point 
spread function (PSF), or impulse response, and tissue is characterized by its reflectivity 
function. With a discrete model for the tissue reflectivity, each scatterer contributes a 
phasor. At any image pixel, with position r e IR 3 , the quantity of interest is the amplitude 
i of the phasor sum, i(r), which has the following form, 

10 ; 

; i(r)=|>Ji(r;r,)e*' (1) 

'=•; where Nr is the number of scatterers in the pixel's resolution cell, and, for each scatterer 
qt is the reflectivity, n = (xi, yi, Zi) is the scatterer position, A(r;n) is the 

15-" position-dependent amplitude of the PSF envelope, and 4h = -2fa& is the position-dependent 
phase where ko is the center wavenumber of the transducer. For simplification, the 
scatterer strength and envelope amplitude can be combined as Ai = qA(r;n) to denote the 
phasor amplitude for scatterer /. 

For any tissue of interest, including the discrete-scatterer surface model, 

20 components of the phasor sum are random. For a given image pixel, the strengths, 
locations, and number of scatterers in the pixel resolution cell are random and produce 
random amplitude, At, and phase, <j>>, in ways that depend on the system PSF. 
Furthermore, for a given surface shape, these interactions change at every pixel, motivating 
the comprehensive model of the present invention. 
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Computation of the mean and variance of the phasor sum amplitude involves 
characterization of the amplitude in terms of probability densities for elements of the ran- 
dom phasor sum (i.e., N, A, and <j>). By assuming that the complex phasor sum is complex 
Gaussian, the sum can be characterized straightforwardly in terms of the system PSF, 
underlying shape and microstructure of the medium. The amplitude mean and variance at 
each pixel can then be computed from the associated parameters for the complex Gaussian 
density. 

B. A Gaussian Approximation to the Complex Sum 

Typically, ./V is assumed to be large enough that the sum is complex 
Gaussian. Formally, this condition holds asymptotically and requires some independence 
assumptions. Practically, the sum can be considered complex Gaussian for N. The sum 
can be written in terms of its real and imaginary components, denoted x(r) and y(r), 
respectively, as 

i(r) = f> eJ0 * W 

= g^cos^+j^Asin^ (3) 
- x(v) + jy(r) (4) 

In general, for a complex random variable, x+jy, the complex Gaussian density is given by 



Pz. y (x,y) = 




where jih and cr x 2 are the mean and variance of the real component, juy and a y 2 are the mean 

E(xy)-E(x)E(y) . . 
and variance for the imaginary component, and r = is the correlation 

(T x <J y 

coefficient. Several assumptions are usually made to simplify this form, but, since some of 
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those assumptions will not always hold, the complete form is retained here. For any 
densities describing N, Ai and ty, the complex sum can be characterized by computing the 
means and variances of the real and imaginary components and the correlation coefficient. 



5 following: 

1 . Independence of individual phasor quantities, Ai _L Aj, (|h 1 <j)j, i * j, 

2. Independence of the number of scatterers from other quantities, 

3. Independence of amplitude and phase for each phasor, Ai _L <j)i , 

4. Identically distributed amplitudes and phases, 

10 5. Symmetric phase distributions that allow rotating the phasor sum to align with the real 
axis. 

For surfaces, assumptions (3) and (5) are too strong to hold at all times, therefore only (1), 
I (2) and (4) are assumed. 

- C. Computing the Complex Gaussian Parameters 
15 The r dependence, indicating the pixel basis for the sum, is dropped from the 

notation for the rest of this paper, although the importance of the implied pixel basis of the 
model cannot be understated. The parameters are given by the following equations: 



The typical assumptions made in the phasor sum analysis involve the 



= E{N)E{A t cos <p t ) (6) 

= £(iY)£(/Min^) (7) 

= E(N)E[A 2 cos 2 <p,) + cr% E 2 {A, cos cj>,) (8) 

= E{N)E{A 2 sm 2 <p l )+a%E 2 (A l sm<j h ) (9) 



\e(N)[E(A 2 cos d> t sin <p z ) - 




20 



All the expressions depend only on the first and second moments of N, A cos 



(fx and Ax sin <jx. Given a shape model, these quantities will generally vary at the pixel level 
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and must be computed from phase, amplitude and number densities derived specifically for 
each pixel. From a computation of these five parameters for every pixel in the image, 
amplitude mean and variance is computed, followed by construction of a data likelihood. 
For surfaces, these results are computed using numerical integration techniques. 
5 The assumptions that are typically made in phasor sum analyses of scattering 

simplify the expressions. First, if amplitude and phase are considered independent for each 
scatterer (At _L (j>0, all expectations involving Ai and cos <|)i or sin (jn can be computed as 
products of expectations, e.g.,E(A i cos <[>. )A i±4 E(A t )£(cos^.) . If, in addition, the phase 
: is assumed to be symmetric about zero, the expected value of sin §i is zero, and the mean, 
10 z- jj, y , of the imaginary component becomes zero as does the correlation coefficient, r. While 
- these assumptions provide convenient simplifications, the results do not hold, in general, 
-. for surfaces. 

D. Computing the Amplitude Mean and Variance 

The amplitude of the complex sum is the quantity of interest in 
15 characterizing clinical images. In the present invention, accuracy in the form of pixel 
densities is weighed against accuracy in describing the differences in pixel densities. Of 
special importance in a model-based approach to shape analysis is how the densities vary 
when the shape changes (i.e., when the likelihood of a different transformation of the shape 
is being assessed). Accordingly, the mean and variance of the amplitude are appropriate 
20 choices for a first approximation. 

For the amplitude of the complex sum, p = ^x 2 +y 2 , the mean, ja P , is 
computed as an expectation over the complex Gaussian density, 
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In general, no closed form solution exists for this integral, thus it is computed for the image 
model numerically using Simpson's rule. 

In contrast to the mean, the variance is quite simple to calculate. Starting 
with the standard expression for the variance, a relationship is stated in terms of the second 
moments of the real and imaginary components and the amplitude mean, 

al = E (p 2 ) - E\p) (13) 
= E (x*) + E (yi) - (14) 

For the image model of the present invention, the amplitude mean and 
variance are the final result. Recall that these quantities are computed from the five 
parameters for the complex Gaussian as in Equations 6 through 10 of this section (i.e., 
section II-B). Those coefficients and the associated amplitude statistics change for each 
pixel. For surfaces, the computations require surface integrals, as described in section III. 
E. The Rayleigh density - an Important Special Case 

For many pixels, the combination of local surface geometry, micro structure 
and system PSF will produce statistics that satisfy the Rayleigh density, parameterized by a 
single constant a, 

p p (p) = -^e-& (15) 

The importance of the Rayleigh density in this work is the need to assign a density function 
to each pixel for construction of a data likelihood. 

The mean of a Rayleigh-distributed random variable is 



ti P = «\/! = \ \J itE{N)E{A 2 l ) (16) 
The variance. <7^, is given by 

a 2 p = f'2-^a 2 = d--^ E{K)E(A? ) (17) 

The simplicity of these relations between the Rayleigh mean and variance and the moments 
of N and A makes the Rayleigh density the choice for the image model of the present 
5 invention. The difficulty is in classifying a given pixel as Rayleigh or non-Rayleigh. In 
this regard, an important quantity is the ratio of the mean to the standard deviation, often 
termed the SNRo, which has a constant value of 1.91 for a Rayleigh density, 

CI SNR 0 = — = 1.91 (18) 

10 * 

Because this ratio is a constant, it is used to classify a pixel as Rayleigh or non-Rayleigh. 
II Specifically, once the amplitude mean and variance are computed (as in the expressions of 
the previous sections), the SNRo is computed, and the result determines whether to 
- approximate a pixel as Rayleigh or non-Rayleigh. Note that a Rayleigh density implies an 
15 SNRo value of 1.91 but that the converse is not necessarily true. 

III. AMPLITUDE MEAN AND VARIANCE FOR SURFACES 
The amplitude mean and variance for each pixel are computed from the 
parameters for the complex Gaussian. The parameters are used to characterize images of 
surfaces. From Equations 6 through 10 of section II-B, the parameters are found in terms 
20 of the local, pixel-based moments for N, the number of scatterers, and the products of 
amplitude, A h and phase, A, cos <fk, and At sin <jx. All of the moment computations 
require surface integrals, where the region of interest is the local intersection of the surface 
and the PSF resolution cell. The resolution cell is defined simply as a region around the 
pixel location in which scatterer contributions are included. Of importance is the 
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intersection surface, Sn, defined as the intersection of the gross tissue surface and the local 
resolution cell. 

A. Computing Moments of N, A and <jt for Surfaces 

The moments of N and A cos (jx require different computational approaches. 
5 Depending on the model used for the number of scatterers, N, the moments for N depend in 
some way on the area of the intersection surface, Sn, 

Arca(Sn) = jj dA (19) 

In Section V, Results, N is modeled as deterministic. With the surface microstructure 
10 parameterized in terms of the scatterer concentration, the moments for N depend only on 
the product of the scatterer concentration and the intersection area. 

The first and second moments for the necessary functions of Ai and (jx require 
computations in addition to the area of intersection. Scatterer amplitude and phase are both 
functions of scatterer position, which is assumed to be uniformly distributed on the surface. 
15 Expected values for functions of that position can, thus, be computed with respect to a 
density for scatterer position. For scatterers uniformly distributed over the surface, the 
density on the scatterer position, r, is simply the reciprocal of the area of intersection, 

pr(r)== r^M for reSn (20) 

The expectation of Ai cos (jx, for example, is then given by a surface integral defined in 
20 terms of position, r, 

E[A,cx>h&,) = jj A, (r) cos o t (r) Pr (r)dA (21) 
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Other moments for products of the amplitude and phase require similar surface integrals 
because the moments are all functions of position. 

Computation of the parameters of the complex Gaussian involves solving 
various forms of the surface integral. The approach selected computes the moments 
5 directly from the triangulated surface. 

B. Numerical Approximations to the Moments 

Numerical approximations to the moments employs approximations to the 
integrals based on the triangles of the surface representation. The area of intersection is 
I approximated as a sum over those triangles in the resolution cell, 

5 Area(S n )= // dA * ^ Area (A,) (23) 

10 ;- • 

The expectation of the amplitude and cosine of the phase is given as follows with the 
« associated numerical approximation, 

K(A,cos&) = jj _ .4,(r)coso,(r) Pr (r)<U (2-1) 

If the amplitude is approximated as constant over the triangle with the value at the triangle 
15 midpoint, Ai(rAi), the computations are simplified further because the integral of the cosine 
of the phase can be calculated analytically, 

E(A, coscM * Area 1 (Sn) A(r A ,)/y co S o,(r)cM 

(26) 

The remaining expectations of functions of At and can be computed similarly with 
20 analytical results for integrals of the trigonometric functions. 

The phase contribution is computed analytically, therefore the accuracy of 
the computations depends on approximating the amplitude as constant over a triangle and 
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approximating the intersection surface by triangles from the original surface. Triangles are 
easily resampled to any desired area, thus the approximations involved (i.e., approximating 
the amplitude as constant over each triangle and approximating the intersection surface by 
triangles) are as accurate as desired without increased computation. 
5 C. Surface Roughness 

The discrete-scatterer model for tissue surface includes a model for surface 
roughness. The surface roughness is modeled as a Gaussian perturbation of each scatterer 
in a direction normal to the surface. The physical perturbation implied by the roughness 
•£ model has not been incorporated into the image model for the sake of simplicity. Note that 
10 f; the discrete scatterer model itself inherently implies a roughness. Adding the physical 
= . perturbation is simple in theory but difficult in computation. If the roughness were 
modeled as uniform instead of Gaussian, the surface integrals of the previous section would 
•i= be volume integrals. The volume integrals would only be significant for the integration of 
j: the cosine of the phase over each triangle and could be calculated analytically. For the 
15 present invention, the additional parameter has little effect since roughness is assumed to be 
small (e.g., a uniform density of widths 0.1 to 0.2 wavelengths). 

VI. METHODS 

In Section V, Results, several images are shown for each of two image 
planes on a vertebra to demonstrate various features of the model. For each image plane, 
20 the following images are displayed: (1) an actual ultrasound image, (2) several images 
simulated for different realizations of the random micro structure, and (3) statistical mean, 
standard deviation and SNRo images. The actual images are registered with the others 
using methods and equipment from image-guided surgery, including optical localization for 
tracking the ultrasound probe and methods for registration of the images and vertebra. The 
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simulated images were generated from the physical model and the statistical images were 
computed in accordance with Section III to characterize all images producible from the 
physical model. 

A. Collection of Actual Images 

A phantom cadaveric L4 lumbar vertebra is shown in Figure 1. The 
phantom was scanned with a CT imaging system to produce an image volume from which 
the vertebra was segmented, allowing the construction of a triangulated surface. A 
rendering of the triangulated surface is shown in Figure 2. Ten aluminum spheres were 
mounted on the outside of the phantom to allow for registration of the physical phantom to 
the CT images. 

Ultrasonic images were acquired using an imaging system from the Tetrad 
Corporation, Englewood, Colorado. A model 6C, 128-element, linear array transducer 
was used. Focus for the transducer is fixed in the elevation dimension and electronic in the 
lateral dimension. Relevant specifications for the transducer include a center frequency of 
6.0 MHz and elevation focus at 33 mm. Based on information from the manufacturer, 
beam width was approximately 1.0 mm in the lateral direction throughout most of the 
image and approximately 3.0 mm in the elevation direction (at a depth of 40 mm). The 
imaging system used, the transducer design and the overall operation of this selection are 
representative of conventional imaging systems. 

The ultrasound probe was modified to allow tracking with an optical 
localization system. The accuracy of using this method to track the position of objects 
identified in the ultrasonic image was measured to be approximately 2 mm. By tracking the 
probe with the same system used to register the phantom and CT images, ultrasonic images 
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were acquired in known relation to the surface model, allowing for a direct comparison 
between simulated and actual images. 
B. System PSF Model 

For the simulations and statistical computations, the system was modeled 

5 using a Gaussian PSF envelope following, A(r;r,) = exp 

The center frequency, fi= 6.0 MHz, was given by the manufacturer. The PSF was 
assumed spatially invariant for this study, and the elevation and lateral PSF widths, oy = 
3 1.5 mm and at = 0.5 mm, were calculated from equations given by the manufacturer. The 
r axial width, oi, was chosen to be 0.2 mm. It was selected from a range of approximately 
10 - 0.15 to 0.4 mm based on visual comparison of actual images (from the Tetrad system) to 
images simulated with the various values for the width. 
r C. Surface Model Implementation - Simulation Procedure 

• = In images of bone surfaces, attenuation is a significant effect that warrants 

H consideration. Therefore, it is assumed that bone is completely attenuating, or occluding. 

15 Before the discrete scatterers are generated, the surface is modified to account for the effect 
of attenuation. A ray-tracing approach is used to determine visibility and is implemented us- 
ing a modified version of a computer graphics algorithm known as Hidden Surface Removal. 
In typical use of this algorithm for rendering a surface to a display, each triangle is projected 
onto the viewing plane, rasterized according to the display grid, and processed under lighting 

20 assumptions to generate an intensity value. During processing, the depth (Z) at each screen 
pixel is stored in a "Z buffer" so that only the closest triangle is displayed. The index of the 
closest triangle is stored instead of the rendering intensity, and, instead of displaying the Z 
Buffer, its contents are used to remove those triangles that were not visible to the transducer. 



{x-Xif {y-yd 2 (z-Zj) 2 



597678.1 4/20/01 



-21- 



An alternative method would be to "clip" those triangles which are partially occluded. This 
alternative is computationally intense. 

After accounting for occlusion, a collection of discrete scatterers are generated 
for the remaining triangles to form the acoustic model of the surface. The distribution of 

5 scatterers are parametrized on the surface according to scatterer concentration 
(scatterers/area) and surface roughness. For each triangle, the number of scatterers is 
calculated as a product of triangle area and scatterer concentration. The in-plane position of 
each scatterer is then generated from a 2D uniform distribution over the triangle. Scatterer 
position, y e IR 3 , is generated from two uniformly-distributed random variables, Aj and X 2 for 

10 triangle vertices, xj, x 2 , x 3 e IR 3 as follows: 

1 . Generate A h A 2 ~ U[0, 1 ] until A] + A 2 <1 

2. y = x 3 + A] (xi - x 3 ) + A 2 (x 2 - x 3 ) 

' where the triangle borders and interior are represented by the combination of the point x 3 and 
- - the vectors x t - x 3 and x 2 - x 3 with A L A 2 e [0, 1] . After the in-plane position is determined for 
15 each scatterer, it is perturbed in the direction normal to the surface to account for surface 
roughness. The perturbation is generated from a Gaussian distributed random variable, with 
roughness characterized by the standard deviation of that perturbation. 

A surface scatterer concentration between 50 and 150 scatterers/mm 2 are used. 
For surface roughness, standard deviations between 0.001 and 0.1 mm are used. Changes in 
20 the concentration and roughness yielded modest changes in the images. The most visible 
change was a decrease in the coherent scattering sites when the roughness was nearly one 
wavelength. In the present examples, scatterer concentration of 64 scatterers/mm2 and 
surface roughness of 0.01 mm were used because texture and intensity produced in the 
simulations were similar to that of the actual images. 
25 D. Computing u, a, and SNRo for Each Pixel 
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The statistics were computed at each pixel in accordance with Section III. 
Surface integrals were computed using the numerical approximations listed in subsection 
III-B. For each pixel, the five parameters of the complex Gaussian were computed as in 
equations 6 through 10 of section II-B, followed by computation of the amplitude mean and 
5 variance as in equations 12 and 14, respectively, of section II-B. Computations are 
performed for each triangle, with the results contributed to the appropriate pixels. 

V. RESULTS 

Statistics computed for all of the image pixels constitute a representation of 
the ensemble of images available from the physical model and, as such, provide a 
10 representation of shape in the images. Statistical images (e.g., m a, and SNRo) provide the 
* basis for a data likelihood but also give a representation of the ultrasonic properties of the 
tissue medium that is fundamentally different from conventional B-mode images. Figures 3 
through 6 represent two different planes on the cadaveric vertebra, as illustrated in Figure 
2. The first image plane shows the transverse process, providing good detail over a small 
15 region. The second plane shows the lamina and articular processes and shows a larger 
region of the vertebra, with a broader range of image features. 

For each image plane, an actual image is shown with sample simulated 
images representative of those possible from the physical model. The statistical mean, 
standard deviation and SNRo images are also shown. The statistical images are computed 
20 directly using the methods of Section III. 

A. A Sagittal Image of the Transverse Process 

The first set of images represent the sagittal plane of the transverse process 
of the cadaveric vertebra marked in the left of Figure 2. The plane corresponds 
approximately to the plane of the actual image in Figure 3. That actual image was acquired 
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by imaging carefully along the transverse process and searching for the brightest (most 
coherent) echo. The sensitivity to angle was such that when imaging by hand, even the 
slightest wobble of the probe changed the image substantially. Consequently, this plane 
represents a challenging surface geometry for the image model since the curvature is 
5 substantial and varied relative to the PSF widths and wavelength. 

In the left column of Figure 3, the actual image is shown above three of the 
infinite set of possible simulated images. Substantial variation exists among the simulated 
images and is representative of the variability expected from the physical model. The 
: actual image appears quite similar to the simulated images in shape, but the actual image 
10 seems to have a wider region of bright echoes along the top of the process with a greater 
*, relative amplitude to the other echoes. Close inspection of the shape in the images shows 
that it is visibly rotated counter-clockwise in the simulated images. Given the sensitivity to 
the angle of insonification, such a change could easily account for many of the differences 
between the images. These differences are unavoidable since the tracking error is on the 
15 order of 2 mm. 

The statistical images show the variation of the statistics across the entire 
image. Dependence on the PSF is evident as the amplitude drops off with distance from the 
surface. A bright region of relatively high amplitude is present in the center of the 
transverse process. From the surface alone, it is difficult to predict that such a region is 
20 produced. This scattering comes from a very small, flat portion of the surface and is 
indicative of the high sensitivity to the surface geometry and angle of insonification. The 
high (greater than the Rayleigh 1.91) value in the SNRo image shows that the region is not 
Ray leigh-distributed . 

B. An image plane along the lamina and articular processes 
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The sagittal plane shown on the right of Figure 2 represents a large region of 
varying surface curvature, featuring both Rayleigh and non-Rayleigh regions as well as 
substantial contributions from the out-of-plane surface in the facet joint (left side of the 
image). Figure 4 shows the actual image and three representative simulated images of the 
5 lamina plane. The lamina and inferior articular process show relatively high-amplitude 
scattering in all of the images, with the lamina producing a higher amplitude more 
consistently. The facet joint is marked by Rayleigh scattering with slowly decreasing 
amplitude toward the top of the image. This slow decrease in amplitude results because the 
41 image plane is nearly parallel to the facet joint surface, with out-of-plane contributions 
10=" decreasing slowly (due to increase in distance from the image plane) over a large axial 
* range. 

= Statistical images for the lamina are shown in Figure 5. Sites of 

I" non-Rayleigh scattering on the lamina and peak of the articular process are evident from the 
-== mean and SNRo images. The greater SNRo along the lamina could be due to an angle of 
15 incidence closer to normal than that of the articular process, although it would be difficult 
to predict by looking at the surface because of the relative flatness of the articular process 
and its wider extent. This ambiguity is another indication of the high sensitivity of the 
interactions between surface and PSF that produce highly variable ultrasonic images. 
Breaks in the non-Rayleigh region along the lamina are due to phase effects and can also be 
20 seen in the simulated images. Similar breaks would be present as artifacts in the 
non-Rayleigh regions if scatterer amplitude and phase had been assumed independent. 
VI. DISCUSSION 
The Rayleigh/non-Rayleigh characterization of surface images does capture 
the basic appearance of surface shape in both the actual and simulated images. Rayleigh 
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and non-Rayleigh scattering occurs with a dependence roughly on the orientation of the 
local surface normal with respect to the axial image dimension. This approach is applicable 
to other tissue surfaces (soft tissue, in particular) and to most conventional imaging 
systems. 

The value of the image model is its basis at the pixel level. That basis is 
especially important for shape because shape must be represented across the entire image. 
The only change at each pixel was the local influence of the shape since the PSF was 
assumed to be spatially invariant and the microstructure homogeneously distributed across 
the surface. In other problems (e.g., phase aberration or attenuation estimation) the PSF 
may change at each pixel. By modeling those changes at each pixel, a likelihood could be 
created that would permit inference of the changes in the PSF or attenuation across the im- 
age. In those problems, a prior model would be expected to provide a substantial boost in 
performance. 

A. Accuracy of computed statistics 

The amplitude statistics employed approximations to the underlying model. 
The computed images shown are accurate in non-Rayleigh distributed regions given the 
high sensitivity to the underlying interactions. An assumption of amplitude phase 
independence, as is typically done in computing density functions for scattering problems, 
would significantly compromise the integrity of the model in the non-Rayleigh regions. 
Visually, the difference of such a choice might be small, but the present invention can be 
directed to the variation of likelihood as assumptions about the underlying structure (e.g., 
the surface shape) are changed. Accordingly, performance in pose estimation for a 
vertebra model is one potential utility of the probabilistic model of the present invention. 
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B. Computational requirements 

The results are generated using algorithms coded in both MATLAB™ and 
C+ + . Computation times were measured for the sagittal-plane lamina images. The 
5 surface region of interest contained approximately 10,000 triangles after resampling by a 
factor of four. The image region was 40x280 pixels. Calculation times listed are for a 
Silicon Graphics Indy workstation. In MATLAB™, computation time for the image 
statistics was approximately 40 minutes per image plane. The same computations 
performed with C+ + code required approximately 8 seconds. The massive speedup is due 

10 to the relative speeds of the two programming environments for calculations involving 
surface triangles and table lookups for approximation of the amplitude mean. The 
computational load is especially important in inference applications since the mean and 
variance must be computed for each image several times (in finite-difference 
approximations of the local gradient of the likelihood) for each iteration of a 

15 derivative-based optimization algorithm. For pose estimation for the vertebra, between 200 
and 300 image calculations were typically required for convergence, when using a 
quasi-Newton BFGS algorithm. 

C. A Data Likelihood 

One objective of the image model was to generate a form capable of 
20 constructing a data likelihood for inference. A data likelihood, p(X\h), characterizes all 
observation data (i.e., all image pixel measurements) with a single probability density 
function, where the data, X = {xi, i = 1, 2, . . . , Nx}, is the set of image pixels, Xi and h 
is some transformation of the surface (e.g., the rotation and translation comprising the 
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surface pose). Assuming that all pixels are independent, the data likelihood reduces to a 
product of probability density functions, pxi(xi), 

p(X\h) = f[p xlh (x i \h) (27) 

From the mean and variance at any pixel, the pixel is represented as either Rayleigh or 
5 Gaussian, depending on the value of SNRo at the pixel. The data likelihood is then a 
product of Rayleigh and Gaussian probability density functions with parameters (oa for xi 
Rayleigh and /m, a t for Xi Gaussian) derived from the system and surface characteristics, 

: ; m>»= ii -v iT n -j==*~ *• " (28) 

x, Rayleigh « ; xfiaussian J7.7ZG "■ 

III where the dependence on the pose, h, is implicit but not shown. 
10; : The likelihood serves the purpose of a cost function, or objective function, 

for estimating the surface transformation. The log likelihood, the logarithm of the 
ri likelihood, preserves the maxima and is simpler to compute in this case because the product 
becomes a sum. 

\n P {X\h)= 2 \*H-4t- S ^ln(2^ ) + ^4^ (29) 

x,^tei g h cc- 2af X]G t^ sian 2 a) 

15 D. Inference 

The representation of shape provided by the image model is illustrated in 
Figure 6. The images in Figure 6 were generated from the Rayleigh/Gaussian-distributed 
image model for the lamina image plane. These images can be compared to the simulated 
images of Figure 4, which were realized directly from the physical model. 

20 The Rayleigh/Gaussian images capture significant information about the 

shape in the image. The relative amplitudes of regions across the image are similar in both 
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sets of images. The differences between the images simulated from the physical model and 
the samples taken from the Rayleigh/Gaussian image distribution are in the image texture. 
Since neighbor interactions are not incorporated in the image model, none of the resulting 
texture should be expected. Inclusion of texture in the image model would be challenging 

5 because of the increase in representation (the joint density over a pixel and its neighbors is 
substantially more complex than the individual density). Application of Markov Random 
Field models could be promising in terms of representing the local effects. While repre- 
sentation of the texture might not be useful in inference of shape, it would be correlated 
with the system characteristics, thus providing a key to inference of those characteristics 

10 (e.g. , the widths of the PSF). 

These and other advantages of the present invention will be apparent to those 
skilled in the art from the foregoing specification. Accordingly, it will be recognized by 
those skilled in the art that changes or modifications may be made to the above-described 
embodiments without departing from the broad inventive concepts of the invention. It 

15 should therefore be understood that this invention is not limited to the particular 
embodiments described herein, but is intended to include all changes and modifications that 
are within the scope and spirit of the invention. 
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